disp(' ')
disp('Descriptive statistics')

% Tabulate the number of ranked school
disp('  Number of ranked schools:')
tabulate(var.nb_sch_ranked_long)

stats.nb_ranked_schools = tabulate(var.nb_sch_ranked_long);

stats.mean_nb_ranked = sum(stats.nb_ranked_schools(:,1).*stats.nb_ranked_schools(:,2)./sum(stats.nb_ranked_schools(:,2),1),1);

% Average % assigned
stats.pct_assigned = mean(mean(reshape(var.stu_assigned_long,I,M),1));
stats.pct_assigned_sd = std(mean(reshape(var.stu_assigned_long,I,M),1));
stats.ex_assigned= numel(var.stu_assigned_long)-sum(sum(var.stu_assigned_long==1));
disp(['  Fraction assigned: ',sprintf('%.5f',stats.pct_assigned),' (',num2str(stats.pct_assigned_sd),')','  (# exceptions: ',num2str(stats.ex_assigned),')']);

% Average % assigned to favorite feasible school
stats.pct_assigned_pref_feas = mean(mean(reshape(var.stu_assigned_pref_feas_long,I,M),1));
stats.pct_assigned_pref_feas_sd = std(mean(reshape(var.stu_assigned_pref_feas_long,I,M),1));
stats.ex_assigned_pref_feas= numel(var.stu_assigned_pref_feas_long)-sum(sum(var.stu_assigned_pref_feas_long==1));
disp(['  Fraction assigned to most preferred feasible school: ',sprintf('%.5f',stats.pct_assigned_pref_feas),'(',sprintf('%.5f',stats.pct_assigned_pref_feas_sd),')',' (# exceptions: ',num2str(stats.ex_assigned_pref_feas),')']);

% Average % students assigned to preferred school
stats.pct_assigned_pref = mean(mean(reshape(var.stu_assigned_pref_long,I,M),1));
disp(['  Fraction assigned to most preferred school: ',sprintf('%.5f',stats.pct_assigned_pref)]);

% Average % students rank top k schools
stats.pct_topk = mean(mean(reshape(var.stu_topk_long,I,M),1));
stats.pct_topk_sd = std(mean(reshape(var.stu_topk_long,I,M),1));
disp(['  Fraction ranked top k_i schools: ',sprintf('%.5f',stats.pct_topk),' (',sprintf('%.5f',stats.pct_topk_sd),')']);